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Abstract 

The long-time behaviour of many dynamical systems may be effectively 
predicted by a low-dimensional model that describes the evolution of a re- 
duced set of variables. We consider the question of how to equip such a 
low-dimensional model with appropriate initial conditions, so that it faith- 
fully reproduces the long-term behaviour of the original high-dimensional 
dynamical system. Our method involves putting the dynamical system 
into normal form, which not only generates the low-dimensional model, 
but also provides the correct initial conditions for the model. We illustrate 
the method with several examples. 
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1 Introduction 



The evolution of many physical systems may be described by ordinary differ- 
ential equations (ODEs) for the "normal modes" of the system. It is often the 
case that most of the modes are strongly damped (these are "stable modes" ) 
while the others ("critical modes") are undamped, or nearly so. In this case, the 
solution of an initial- value problem rapidly approaches a low-dimensional centre 
manifold (denoted by M), which may be parameterised by the amplitudes of 
the critical modes. The subsequent evolution of the system on J\A is slower than 
during the initial rapid approach to M.. 

The centre manifold is an important theoretical and applicable tool for sev- 
eral reasons. Firstly, we arc often concerned with the stability of a fixed point 
to small disturbances. The stability and bifurcations of the fixed point may 
be completely determined by analysing just the low-dimensional evolution re- 
stricted to M. H Q . Rational approximations to M. may be computed, that is, 
there exist constructive schemes for approximately eliminating the stable vari- 
ables || [17| [18], [l9], go], |2^, |25|, |2(| . Such schemes involve approximating M by 
(the first few terms in) a power series in the amplitudes of the critical modes. 
Secondly, if the fixed point is stable then to any solution Pit) of the origi- 
nal problem there corresponds a solution Q(t) on M. which P(t) approaches 
exponentially quickly. The long-time behaviour of the full problem is there- 
fore determined by the low-dimensional (and hence tractable) centre-manifold 
model. 

The subject of this paper is the relationship between the full solution Pit) 
and its approximation Q(t) on M. Specifically, we consider how to derive the 



initial value Q(Q) from P(0). For although the theory of centre manifolds guar- 
antees the existence of Q(t), it gives no general constructive means of determin- 
ing Q(0) (a method for small -P(O) — Q(0)| was given by Roberts |24]). Once 
Q(0) is known, the solution Q(t) at later times is computed by integration of a 
low-dimensional system of ODEs. 

Computing Q(0) by a perturbation expansion in a small parameter related 
to the time-scale of the approach to the centre manifold has been discussed by 
van Kampen Q, see also the references therein]. However, his method requires 
the solution of a succession of ordinary differential initial-value problems, which 
is often impractical. In Sections 2 and 3 we develop a simpler procedure for 
calculating Q(0), by algebraic manipulations alone: no ODEs need be solved, 
as we illustrate in a simple example in Section 4, where we compare our method 
with van Kampen's. The heart of our method is putting the evolution equations 
into normal form [g, [7] . A well-known by-product of the reduction to normal 
form is the centre manifold fa] . The novelty of our work lies in the observation 
that appropriate initial conditions for the centre-manifold model follow naturally 
from the normal-form reduction. 

We further illustrate our ideas in Section 5 with Taylor's model of shear 
dispersion in a channel, where both the original system and the model are partial 
differential equations. Nevertheless, the model, which involves just one spatial 
dimension, is considerably simpler than the original problem, which involves 
two space dimensions. However, at present there is little theory to support 
an infinite-dimensional centre manifold, so the treatment is necessarily formal. 
As a further extension of our approach, we discuss in an appendix the use of 
normal form transformations in constructing initial conditions for dynamical 
models that involve both unstable and stable modes. 



2 A system with centre and stable modes 

We consider in this section a nonlinear system of the form 

x = Ax + M(x, y) xefi m (1) 

y = By + N(x,y) y £ M n , (2) 

where the eigenvalues a& of the m x m matrix A satisfy 5R(a/c) = 0, and the 
eigenvalues (3k of the n x n matrix B satisfy 5R(/3fc) < —ft < 0. The functions 
M and N are strictly nonlinear smooth functions, of O (|(x, y)| 2 ). 

Under these conditions centre manifold theory may be applied Q to deduce 
the following results. 

CM1 There exists an m-dimensional centre manifold M. of the form y = h(x), 
with h(x) = O (|x| 2 ). On M the m-dimensional dynamics of x are de- 
scribed by ([I]) restricted to A4, namely 

x = Ax + M(x,h(x)). (3) 



CM2 Corresponding to each solution P(i) = (x(i),y(i)) of p H) there is a 
solution x = u(£) of (||) , such that 

\P{t) - Q(t)\ = O (e- pt ) asi^cx), (4) 

where Q(t) = (u(t),h(u(*))). 

A sufficient condition for this property to hold is that the origin should 
be stable. This condition is not, however, necessary, and we assume that 
CM2 holds regardless of the stability of the origin. 

CM3 If an approximation 0(x) to h(x) satisfies 

50 + N(x, 0) - ^ L4x + M(x, 0)1 =0(|x| p ) as Ixl -> (5) 
ax 

thcnh = 0(x) + O(|x| p ). 

Solutions of (|l]-||) therefore are exponentially attracted to M. (at least for suffi- 
ciently small initial values of x and y) , and at large times the essential dynamical 
behaviour of (|l]-[|) is captured by (Q). This is an important result for practical 
purposes because (g) is of low dimension, m, compared with (flj-g), which is of 
dimension m + n. Indeed, centre manifold theory also applies to problems with 
infinitely many stable modes H, for example, the Kuramoto-Sivashinsky equa- 
tion pi. Therefore, when interested only in the long-term behaviour of (ffl-0) we 
may as well compute solutions Q(t) of (ph rather than solutions P(t) of (]l]-g). 
Furthermore, an important property of (1TI--0) is that if |/3j| S> 1 for some of the 
eigenvalues (3j then the system is stiff, which makes reliable numerical computa- 
tions expensive. The model (pi) is free of such undesirable stiffness (although it 
may have the different problem of rapid oscillations), which facilitates numerical 
computations, particularly over long times. 

In general h(x) cannot be computed exactly. But by a straightforward al- 
gorithm involving algebraic manipulations based on (H), a power series in x 
may be developed for h(x) a, Elf. In the following subsections we describe a 
similarly straightforward algorithm, based on a normal form transformation, for 
determining the correct initial value Q(0) from the initial value P(0), so that the 
exponential approach (0) is achieved. Only with a correct initial condition can 
long-term quantitative predictions be made by the low-dimensional model (||). 

2.1 Normal form and centre manifold 

Our method for constructing appropriate initial conditions relies on the struc- 
ture of the normal form transformation. As shown by Elphick et al. 0, for 
example, it is always possible to find a nonlinear change of coordinates, to vari- 
ables x € lR m and ?7 e iR", of the form 

x = x + F{x,v)v (6) 

y = v + s(x,v) (7) 



such that the dynamical evolution of (UH3) may be written in normal form 

X = A x + a(x) (8) 

r) = B V + B(x,ti)v- (9) 

In these equations: F is an m x m matrix, B is an n x n matrix, and both 
are at least linear functions of their arguments. The functions g € M n and 
a G M m are strictly nonlinear. The nonlinear terms that remain in the normal 
form equations (pH|) arise from resonances between eigenvalues of the linear 
operators |g|. 

Note from (0) that the centre manifold 77 = is invariant under the evolution 
of the system. For small x and r\ the linear term Br\ dominates the nonlinear 
terms in (m, and so 77 — > exponentially quickly. In terms of the original 
variables, M is given by 

x = x, y = g(x,o), 

or y = h(x), where h(x) = g(x, 0). Equations ([|) and (g) describe the same 
long-term dynamics. 

2.2 Initial conditions 

We now show that the normal form transformation gives as a by-product an 
appropriate initial condition Q(0), so that we can make long-term predictions 
with (|). 

First we observe that if the initial value Xo is fixed, all solutions of the 
transformed equations (g |9j) from initial conditions (xo, ?7o) have the same long- 
term dynamics, regardless of the value of 770 (provided it is small enough to 
guarantee approach to M). This is because the evolution equation (||) for x. 1S 
independent of r\. Therefore (x j 0) is the initial condition on M. for a solution 
that soon evolves identically to the solution from the initial condition (xo;?7o) 
off M. 

Now it is clear how to derive an initial condition Q(0) for the model (||) from 
the initial condition (x(0),y(0)) = (xo,yo). 

1. Use the coordinate transformation (jq-wj) to determine the corresponding 
initial condition {xo, 770) for the normal form equations. 

2. As noted above, the corresponding initial condition on M. is (x(0), 77(0)) = 
(Xo,0). 

3. By the coordinate transformation (Q), the initial condition for the model (||) 
is then x(0) = xo; that is, Q(0) = (xo,g(xo, 0)). 

Since the change of variables is determined in powers of the new variables x 
and 77, the evolution equations of the normal form, (@) and (H), are equivalent 



to (|l|-§) only up to terms smaller than any power of |(%, 77)! as | (x, ^7) | — > 0. 
The differences will become significant if our initial point is close to the bound- 
ary between the basins of attraction of different attractors on A4 (particularly 
this will be a problem if the attractors have fractal basin boundaries) or if the 
solution on M. has a positive Lyapunov exponent. In these cases our method 
may, like almost any other approximation, fail reliably to predict the long-term 
behaviour of the system. An analysis of this limitation is beyond the scope of 
this article. 

Provided the initial value of r\ is not too large, all initial conditions for (|8|- 
|9|) in a given plane of constant \ evolve in essentially the same way at large 
times. In the original (x, y)-space, the planes II(xo) = {(X; v) I X — Xo} appear 
as curved manifolds S: 

s (xo) = {(x, y) I x = xo + F(xo, v)v, y = V + g(xo, v)} . 
where \q is fixed on each manifold, and r\ parameterises it. 



The manifolds S(xo), which Roberts 24 has termed "isochronic manifolds", 
are a generalisation of the concept of "isochrons" introduced by Winfree |U[ 
in the context of resetting biological clocks. Winfree considered a nonlinear 
oscillator with a stable limit cycle T as a model for a biological clock. Each point 
on r can be assigned a phase <fi which increases uniformly in time, <f> = u. If a 
solution on the limit cycle is perturbed away from T then it will relax back to T— 
the biological clock will reset itself — but with a slightly altered phase. Winfree 
proposed the concept of isochrons to describe how different disturbances induce 
different changes to the phase of the solution. To each point Q on T, with phase 
<j> say, he associated a surface £(</>), transverse to T, passing through Q. This 
isochron £(</>) consisted of all points P which, after their transient approach to 
F, have the same phase as Q has at the same time. Later, Guckenheimer [O 
proved the existence of isochrons for a nonlinear oscillator, using the invariant 
manifold theorem. 

Our method for computing appropriate initial conditions for a centre-manifold 
model, is founded on the fact that projection along the planes II is simpler than 
nonlinear projection along the curved surfaces £. 



2.3 An example 

Consider the system 



X = XQ.-X-0Y) (10) 

Y = Y(l-Y-aX), (11) 



which models a range of phenomena from population dynamics to competing 
modes near a multiple bifurcation point. To demonstrate our method we exam- 
ine the degenerate case of (3 = 1, and make the linear change of variables 

X = x, Y = 1 + y — ax 



in order to bring (^-0) into the form (§-[2|), with 

x = -(l-a)x 2 -xy (12) 

V = -y - a(l - a)x 2 - y 2 . (13) 

We now apply the first nonlinear stage of a normal form transformation to 
remove unnecessary quadratic terms from (O-O), by setting 



x = x + ®X 2 + bxv + crj 1 (14) 

y = V + dx 2 +e\V + fv 2 , (15) 

with constants a-f to be chosen. Under this change of variable (]12|-|l3|) become 

x-bxn-zwf = -(! -«)x 2 -xv + o(\(x,v)\ 3 ) 

il-exr\-2jif = - V -[d + a(l-a)}x 2 -exV-[f + Mv 2 + 0(\(x,v)\ 3 )- 

To simplify the equations that govern x an d V we choose 6 = 1, c=0, d = 
—a(l — a), / = 1. The coefficients a and e remain at our disposal and we 
choose to set them to zero. The evolution equations for x an d f] are then, to 
quadratic order, 

X = -(l-a) X 2 (16) 

17 = -r,. (17) 



Suppose now that we are given the initial condition P(0) = (xo, j/o) f° r (12 
|13| ) near (xo,J/o) = (0,0). Then by (|14|-jl5|) the corresponding values for x and 
r\ are 

Xo = %o -Xoyo + O(\(x ,y )\ 3 ) 

»?o = J/o + a(l -a)x„ - j/q + (K^o^o)! 3 ) ■ 



According to (flTn 7/ — > exponentially as t — > oo, and x evolves independently of 



77, so the long-term evolution in (16-17J) is the same as from the initial condition 
(X0)0). This corresponds to the initial condition Q(0) — (xq,j/o) f° r (Ftt^D? 
where from (nJ-OS) 



Xo = xo - x y 

„2 



Vo = -a(l-a)x . 

3 Slow-manifold models 

Suppose now that the rapidly decaying modes of (0-0) have been eliminated, 
and that the dynamics have been reduced to (0) on the centre manifold. Some 
eigenvalues of A represent slowly evolving modes, others fast oscillations. In 



this section we discuss models in which fast oscillations are eliminated, as the 
rapid transients were eliminated in the previous section. This is the essence of 
many physical approximations, and is a useful tool for numerical integrations. 
This procedure finds application in fluid mechanics, where high-frequency sound 
waves are ignored under the incompressible approximation [JI4], §9]; in beam 
theory, where fast ringing modes are ignored |26| ; in meteorology, where weather 
data must be "initialised" before numerical forecasts can be computed, in order 
to remove spurious, relatively high-frequency gravity wave oscillations || |l5|, [l6| . 
Suppose the linear spectrum of A consists of a zero eigenvalue repeated m 
times, together with n/2 complex-conjugate pairs of purely imaginary eigenval- 
ues ±ia>fc. Then with some renaming of variables we may split (|3|) into slow and 
fast components, 

x = Ax + M(x,y), y = By + N(x,y), (18) 

where x € M m , y 6 M n , A has m eigenvalues that are precisely zero, and B 
has purely imaginary eigenvalues, {j3\ , . . . , /3„) with /?2fe = —fak-i- Sijbrand 
[ pT| shows that, since there is no resonant forcing by x of the fast oscillations 
in y, this system has a sub-centre manifold M.s of the form y = h(x): a "slow 
manifold" that may be derived for this system by the same formal means as 
the centre manifold described earlier. However, the slow manifold does not in 
general attract neighbouring solutions, but instead acts as a slowly evolving 
"centre" for their rapid oscillations. 

For a slow-manifold model the normal form for x cannot in general be made 
independent of 77 — it takes the form 

X = A X + a{x,v), (19) 

with a(x,7?) = 0(|?7| 2 ) as \r/\ — > 0. Quadratic terms in 77 arise from resonant 
forcing of the slow modes by nonlinear combinations of two conjugate modes 
±iuJk- The evolution equation for the fast variables 77 may still be written in the 
form (||). The subcentre manifold devoid of the fast oscillations characterised 
by the imaginary eigenvalues of B is given by 77 = 0; the evolution of the slow 
variables x is given by (|l^) restricted to 77 = 0. 

3.1 Initial conditions 



Integrating <fiM numerically can be costly if the rapid oscillations (potentially 
of large amplitude) are present. Unfortunately, since (M) cannot in general be 
made independent of 77, solutions on and off the slow manifold (where 77 = 
and 77 7^ 0, respectively) evolve in essentially different ways. It is therefore 
impossible in general to project an initial condition (xo,yo) to Ms so as to 
maintain the same essential dynamical behaviour for all time. However, it is 
possible to match the behaviours for a time of o(|t7|~ 2 ) in the presence of fast 



oscillations of small amplitude \r]\. This approximation, where only those terms 
up to linear order in 77 are kept, is described by Roberts |24|]. 

A more favourable situation occurs when the system (|18|) is symmetric. If 
M(— x, — y) = — M(x, y), and similarly for N, then a fast mode cannot resonate 
with its conjugate alone to force the slow modes. Instead a combination of at 
least three modes of at least two different frequencies is needed to force the slow 
modes. Hence in this case a(x,i]) — O (\rj\ q ) for some q > 3 and the effects of 
any fast oscillations on the long-term evolution are much weaker. 

3.2 A simple example 

Consider the system 

x = yl, 2/i = -2/2, 2/2 = 2/1, 



which is of the form (18). The y- variables execute fast oscillations with period 
2n, while the variable x is constant according to the linearised equations, but 
is forced quadratically by the rapid oscillations. Here the normal form trans- 
formation can be accomplished exactly. We set \ = x + 2/12/2/2, 771 = y%, and 
7/2 =2/2, so that the normal form equations become 

x = (vi + vl) A »7i = -»72) m = vi, 

or 

X = r 2 /2, f = 0, = 1, (20) 

where we have written 771 = rcosO and 772 = rsin#. Thus r is the amplitude of 
the oscillations, and 6 is their phase. From (EOJ), the invariant slow manifold, 
Ais, is given simply by r\ = (i.e. r = 0); the evolution of a solution Q(t) 
on Ais is trivial: % = 0. So solutions Q(i) have \ constant and 771 = 7/2 = 0. 
However, for solutions P off Ais, r is a non-zero constant, and so x drifts at a 
(non-zero) constant rate: x — r 2 /2. It is not possible to match solutions P(t) 
and Q(t) for large times by any choice of initial condition (3(0). 

4 A comparison of our method with others 

In order to show how the method we have described for determining P(0)—Q(Q) 
compares with other methods, we describe now three alternative approaches for 
a simple problem. (The difference P(0) — Q(0) is sometimes termed "initial 
slip" prt by analogy with the boundary slip of an inviscid fluid.) These are: in 
the first instance to solve it exactly; secondly to apply our normal-form method 
described above, and finally to use explicit perturbation expansions for the 
solution (rather than for the governing equations, as is the case for the normal 
form method) to derive the same results. In general, for nonlinear problems, the 
first approach is not available to us, and the algebraic details of the third quickly 



become overwhelming. A further method for computing the "initial slip" , which 
we shall not discuss here, has been derived by Geigenmuller et al. M, who apply 
a systematic perturbation procedure to linear systems of a certain form. 

The example involves a simple chemical reaction [H4L §2] which in scaled 
variables that represent two concentrations is written 

(21) 
(22) 

In the first equation x evolves slowly, because we consider e to be small. These 
equations are linear and so the system may be solved exactly, in particular the 
behaviour of the system at large times is determined exactly. 

4.1 Exact solution 

Equations (pl|-p2|) may be rewritten as 



x = -e(x-y) 

y = kx — y. 



d 


X 


= M 


X 


dt 


[ y _ 




[ y \ 



where M 



— € € 

k -1 



(23) 



The eigenvalues Ao and A s of the matrix M are 

Ao, A s = 1 (-(1 + e) ± V(l-e) 2 +4efc) , 

As e — > 0, the eigenvalues take the form Ao ~ e(k — 1), A s ~ ■ 
to slow evolution and rapid decay, respectively. 
The solution of (§|) is 



-1, corresponding 



x(t) = ae A °* + be 



„A a * 



y(t) 



ka 



„Aot 



kb 



^x s t 



1 + A 1 + A s 

where the constants a and b are found from the initial values: 



(1 + A )(1 + A S ) 



-y + kx /(l + \ s ) 
y - kx /(l + A s ) 



fc(A - A.) 

Since 1 3> e(fc — 1), the A s -mode decays to zero much more rapidly than the 
Ao-mode evolves, so that after some time the solution is given, to within expo- 
nentially small terms, by 



x(t)~ae Aot , y{t) 



ka 



,A t 



k 



(24) 



1 + A 1 + A 

This solution is equivalent to the solution of the one-dimensional model equation 

X = X X (25) 

with initial condition 

(1 + A )(1 + A S ) ( , kx 



where x = X and y = kX/(l + A ). 



-J/o 



1 + A, 



(26) 



10 



4.2 Normal form calculation 

The normal form for (pi|-|22|) is obtained by making the substitution 

x = X + fv, V = V + 9X, (27) 

where /(e) and g(e) are constants. Wycoff & Balasz [B2( have considered lin- 



ear systems of a form that includes (01-22), and have derived a substitution 



of the form (E7f) by considering a multiple-time-scale perturbation expansion. 
Our results for this problem agree with theirs for the more general case. Equa- 
tions (^I|-f2"2|) fail to satisfy the conditions of (|l] ||) in two ways: firstly, there 
is no mode with precisely zero growth rate; secondly, the system has not been 
decomposed into critical and stable subspaces. The first point is dealt with by 
using a standard trick of bifurcation theory H , namely considering the param- 
eter e as an extra dynamical variable with the trivial evolution equation 

e = 0. (28) 



Then the term — e(x — y) is nonlinear in the extended system (|2l| , |22| , |28| ), and so 
the linearised evolution of x is just x = 0. The second point is unimportant in 
practice, and we diagonalise the system as we go. 

Substituting (E7fl into the governing equations (EII-E4) we find 



x + fv = 


-ex - efrj + en + eg\ 


V + 9X = 


-V - 9X + k X + kfn. 



Separating out terms proportional to x an d f] we find that the slow and fast 
variables evolve according to 

X = -e(l-g)x (29) 

f, = -(l-kf)n, (30) 

where / and g must satisfy quadratic equations whose roots are 

/ = i (l - e - v / (l-e) 2 + 4efc) ~ -e as e - 
1 



g = — [ e - 1 + v / (l-e) 2 +4efcj ~ k as e -> 0. 

The fact that g is not small as e — > reflects the fact that the original system 
is not diagonalised. 

The centre manifold may now be read off from the normal form (|27|) by 
setting 77 = 0. That is, 

x = X, V = 9X- (31) 

Since g(l + Ao) = k, this expression for the centre manifold agrees with (p3) 
derived from the exact solution (as it should). Since — e(l — g) = Aq, the 



11 



evolution of the slow variable \ m (P9|) i s precisely that derived from the exact 
solution, where we denoted the slow variable by X in (|25|). 



Now we follow the procedure described in Section 2.2 to determine the ap- 
propriate initial condition Q(0) on the centre manifold which gives rise to the 
same long-term dynamics as P(0) = (xo,ya). Using ( p7| ) we find that -P(O) is 
written in terms of \ an d rj as 

(Xo, Vo) = -, ^(^o - fyo, -gx Q + y ). 

The expression for the corresponding Q(0) on M. is (xo>0)- Transforming this 
point back to the original variables yields the expression 

{%o,yo) = — -j- ( x o - fyo,g[xo - fyo]) ■ (32) 



A little algebra shows that ( |32| ) agrees with the exact solution (|2 

4.3 Perturbation expansion of the solution 

Finally, we treat the problem with a techniquebased upon a perturbation ex- 
pansion for the solution during its rapid approach to the centre manifold. Such 



a method is described by van Kampen 1 14 for problems that may be nonlinear, 



and we follow his argument here. The essential idea is that during the rapid 
approach to the centre manifold, the fast variables change a great deal, while 
the slow variables change little. 

We begin by expanding the slow variable x(t) as a power series in e, 

x{t) = x°(t) + ex l {i) + e 2 x 2 {t) + ■■■, (33) 



and then we substitute this expansion into (|2l|-|22|). Note the important dif- 
ference between this perturbation expansion of the solution and the expansion 
of the governing equations for the normal form calculation. Approximations 
to the fast variable y(t) are constructed from (|22|), with x replaced by a finite 
truncation of the series ([53|) . 
At O(e ) in @, ±° - 0, so 



.r 



= constant = xq. (34) 



To leading order then, during the approach of the solution to Ai, i does not 
change. Now we compute an approximation to y(t) during the approach to M. 
using © and (H), 

y = -y + kx°, 

which may be solved to give y(t) = yoe~ l + kxo(l — e - *). Now we usethis 
solution for y(t) to compute x\, by considering (Elf) at O (e 2 ): 

x 1 = -x°+y(t), 
12 



whose solution is 

x\t) = -x t + y (l - e _t ) + ky (t + e - * - 1). 

Now we re-assemble the solution x(t) to O (e) and find 

x(t) ~ {x + e(y Q - kx )} + et {-(1 - k)x Q } + ee~* {-y + kx } ■ 

The first two terms are the start of a Taylor expansion of the solution in powers 
of £, while the last term is an exponentially decaying transient. The initial value 
for the Taylor-expansion component, xq + e(yo — kxo) , is just the initial value for 
the solution Q(t) on M.. This result agrees with the previous two calculations 
of the initial condition Q(0), as we see by expanding ( p4 ) to give 

* xp - fy x + ey Q . 2 , 

x a = 1 _ . ^7^ = ^o + e(?yo - ^o) + C (e ) . 

To continue the expansion requires the solution of a succession of ODEs: in 
general, even if possible, the computations rapidly become exhausting. 

5 An infinite-dimensional example: shear dis- 
persion 

We now apply the technique we have described for determining initial conditions 
to a generalisation of Taylor's model p8f for shear dispersion. This model 
idealises the spreading of some contaminant released into a river. The river is 
modelled as a channel, and the water is supposed to flow parallel to the banks, 
with the flow being fastest closest to the centre of the channel. The variation 
in water speed with distance from the bank increases contaminant gradients, 
while molecular diffusion tends to smooth gradients out. A balance between 
these competing physical mechanisms occurs at large times after the release of 
the contaminant; eventually the contaminant is approximately evenly spread 
in any cross-section of the river. The cross-sectionally averaged concentration 
varies slowly in the downstream direction. The slow evolution of the averaged 
concentration obeys a simplified model equation which was derived through the 
techniques of centre manifold theory by Mercer & Roberts |17], |l8| . 

The following analysis can be made rigorous in Fourier space p7| , but this 
rigour is generally unavailable in other problems where the centre manifold itself 
is of infinite dimension. In the absence of rigorous theory, the best we can do is 
to apply our techniques formally. 

5.1 The model and its linear behaviour 

Consider the dispersion of a contaminant, with concentration c(x,y,t), in a 
channel < y < 1, — oo < x < oo. If longitudinal spatial variations (that is, 
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variations in x) occur on a large scale than c satisfies the partial differential 
equation (PDE) 

dc , . dc d 2 c ,„ . 

m = - u ^T x + W (35) 

subject to the boundary conditions that no contaminant escapes through the 
walls of the channel, 

-^=0 at y = 0,1. (36) 

dy 

The variables have been made dimensionless, and u(y) is the velocity in the 
2-direction. 

The PDE ( p5|) has two space dimensions; the "centre manifold" we derive 
by performing a normal form transformation involves just one space dimension. 
We concentrate on the structure of the normal form transformation; the details 
have been discussed elsewhere JL7], fL8|. 

The problem (|35|-|36|) is linear in the concentration. The concentration 
evolves so that after some time the spatial variations in x are small, and so 
we treat all ^-derivatives as "nonlinear" terms in the same way that a small pa- 
rameter was treated earlier. With this interpretation, the "linearised" dynamics 
are given by 

dc d 2 c 

dt dy 2 ' 

which physically describes cross-channel diffusion. The eigenmodes e ant c n (y) 
are of the form 

c n = cos nny (n = 0, 1, ...), 

with corresponding eigenvalues a n — —n 2 ir 2 . According to this "linear" picture, 
all modes decay except cq. The system is therefore analogous to (0-||), with one 
critical mode but infinitely many damped modes. Note that the critical space 
is spanned by the single mode, Co, and that all other eigenmodes have zero y- 
average. For a given concentration field c we can determine the component in 
the critical space by taking a y-average to give c; then the component in the 
stable space is c' = c — c. 

5.2 Normal form 

Here we derive a "normal form" for (pBJ). The normal form transformation 
consists of choosing new variables x an d V m terms of which the governing 
PDE (|35| ) is simplified. The general form of this transformation is 

c = x + 9i(x,v), c' = r) + g 2 {x>v)- 

Since the governing equation is linear in c then so are <?i and <?2- They 
are "nonlinear" in the sense that they vanish when x-derivatives are ignored. 
Using experience from the previous sections, we choose c to reduce to x when 
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r\ = 0, and c' to reduce to r\ when \ = 0. Furthermore, we suppose that x 
has no component in the stable space (x = \) and that r\ has no component 
in the centre space (77 = 0). Such considerations imply that the "nonlinear" 
terms in the normal form transformation must be of the form 31 (Xi 7 ?) = Sr/, 
92 ix-iV) = h-Xi where h and 9 are linear operators such that 9rj = 9r\ and 
h\ = 0. For definiteness we define 9 so that 6>/ = for all /. We have therefore 
decomposed c(x, y, t) as follows: 

c = x + hx + W + 0v- (37) 

Because the original problem is linear in c we may deal with the x~ an d 
the 77-components separately. After substituting (|37]) into ( |3q ) we consider first 
those terms involving x- We find that 

{l+h)d l = - u h {l+h)x] + w hx - (38) 

This equation serves to define both the operator x and the derivative dij/dt. 
The average of this equation over y gives an expression for dx/dt: 

This is the evolution equation for x on the centre manifold, and represents a 
simplification of the original problem (|3|). Not only is ( |39| ) of lower dimen- 
sion than (p5|), because there is no y-dependence in x(x,i), but there are only 
"nonlinear" terms in (B9h — that is, every term on the right-hand side of (JM) is 
differentiated at least once with respect to x. 

Substituting (B9h back into (B8), we obtain an equation that governs the 
operator h: 

— - -ufl h)— — 
dx dx dy 2 

In order to solve this equation, we expand h in powers of d x , 

00 

j=i 

substitute this expression into pQ), and equate powers of d x . Retaining only 
the first term in the expansion of h, we find ( p9| ) to be approximately 

dx -dx d 2 x 

at ox ox z 

The concentration x is advected with the mean velocity, u, and diffuses with 
diffusion coefficient —ua\. 
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(1 + h)u{ i + /,)_ = _ tt(1 + ft )_ + h . (40) 



A similar argument for the terms involving 77 yields 

(1 + e) |._„| [(1 + ,), ] + 0. 

Averaging with respect to y gives 

Subtracting one equation from the other gives an evolution equation for 77: 



and so 6 satisfies 



® n\ i d 2 r]\ d 



This equation may be solved by expanding 9 in powers of d x . Note that the 
manifold r\ = is invariant (from (pd|)) and the evolution equation (41) for 77 



contains both "linear" and "nonlinear" terms. Since rj = then the eigenvalues 
of the "linear" operator are Ai, A2,. . . , which are all negative. Consequently, 
?y — > exponentially quickly. 

If we compute the terms in the expansions for h and 6 we find that Qr\ = —hr), 
so the normal form transformation (p?|) puts c into the form 



c = x + hx + T) — hr]. (42) 

Note that when we apply the operator 1 + h to (|4|) and average with respect 
to y we obtain 

(1 + h)c = JTTWx- (43) 

This enables us to determine \ from a given c. To then determine the corre- 
sponding value for 77, we use 

77 = c — c — ft-x- (44) 

5.3 An appropriate initial condition 

To derive initial conditions for x so that its long-term evolution in (|3^) agrees 
asymptotically with that of c in (pq ) from a given initial condition c= Co(x,y) we 
first apply (f43J-[44|) to obtain the values (xo, Wo) that correspond to cq- Then we 
note that the solution of (|39|~[4l|) from P(0) = (xo, Wo) exponentially approaches 
the solution from Q(Q) — (xo,0). Transforming Q(0) back to the physical 
variable c, using (Eq), we obtain a modified initial condition, cj$, for (B5|) that 
lies on the centre manifold and has the same long-time behaviour as from Co. 
Of course, if we want to integrate the model equation (p9) instead of (pq) itself, 
then all we need is the value xo- 
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5.4 Comparison with systematic adiabatic elimination 

In order to emphasise the relative simplicity of our method for determining 
"initial slip" by purely algebraic manipulations, we now describe an alternative 
method (Haake and Lewenstein [l^] ): that of systematic adiabatic elimination. 
(See also a treatment of ( pq ) in a different context by Titulaer p9[.) The 
governing equation (pq) for the concentration c(x,y,t) is written in the form 

dc 

— = {L + Li)c=Lc, (45) 

where for the shear dispersion problem 

Formally, the solution may be written in terms of the initial concentration dis- 
tribution as c(x, y, t) = e Lt c(x, y, 0). 

An evolution equation is sought for a "reduced" variable C(x, t) — c. This 
equation will be of the form 

at - £C ' 

for which the formal solution, subject to the initial condition C(x,0) = C c g(x) 
is C(x, t) — e et C c g(x). Our goal is to calculate the operator £, and to determine 
C e ff (x) so that the reduced description of the dynamics matches the full solution, 
that is, 

e Lt c{x,y 7 Q) ~ e H C eS (x) 



for large times. In the notation of the previous subsections £ = —u(l + h)d x 
and Ceff(x) = x(a;,0). 

The computation of these quantities by systematic adiabatic elimination 
proceeds as follows. First we use (|4^) and the definition of C to write 



C(t) = c(t) = c(0) + / Lc(r)dT = c(0) + / Le LT c{0)dr, (46) 

where we have suppressed the dependence of c and C on x and y for notational 
brevity. Now we assume that d Xl and so L\, is small, and we carry out a small-i, 
small- Li expansion for C . The first step is to approximate exp(Lt) by exp(io^)- 
Then we note the decomposition 

e Lo *c(0) = c(0) + [e Lo *c(0) - c(0)] . 

The term in brackets is just 



L^-e Lat c(0), 
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where the "inverse" of Lq has been taken to have zero y-mean. Therefore, 
collecting these results and substituting them into (E6h, we find 



ft f°° f) 

C(t) « c(0)+/ L lC (0)dT+ L 1 Lo 1 —e L ° T c(0)dT 
Jo Jo ut 



« c(0)+iLic(0)-LiL ( 7 1 c(0) + ---. 

Here we have extended the range of one of the integrals to infinity, and incur 
exponentially small errors as a result. This expression represents the Taylor 
series expansion of C(t) for small t, with initial condition 



C eff (0) = c(0) - LiLq x c(0). 

= ^) + ^u(y)Lo 1 c(0) 
d 



= c(0) + —c(0)L^u(y). 

This result agrees with the previous calculation (where we use the result that 
ai = i^Vy)@,|l|). 

To solve for even this leading-order contribution to the initial slip (there are 
also contributions from all higher powers of d x ) we had to evaluate an integral 
of the exponential of some operator. For shear dispersion, the integral can be 
evaluated explicitly, although at higher orders and in truly nonlinear problems 
the exact computation of such an integral will not be possible. By contrast, 
the method we have proposed for computing initial slip by purely algebraic 
manipulations is much more straightforward. 

6 Conclusion 

We have described a method that uses a normal form transformation to de- 
termine appropriate initial conditions for low-dimensional models of evolving 
systems. The nonlinear coordinate transformation allows one to compute the 
initial slip by solving a succession of algebraic problems, rather than the differ- 
ential problems associated with other methods. The process has been illustrated 
with several examples. Using a normal form transformation is simpler than other 
methods, and also shows, for example, how nonlinear resonances make it im- 
possible, in general, to find good initial conditions for slow sub-centre manifold 
models such as the quasi-geostrophic approximation or beam theory. 

Generally, the calculation of the initial slip will be carried out only approxi- 
mately (usually by resorting to the truncation of a power series), so the solutions 
of the model and the full system will differ to some order, after a period of rapid 
exponential approach. Since the change in variables is made as a power series in 
the new variables, the approximate normal form equations differ from their ex- 
act counterparts by, say, terms of order ((%, rj)\ as \{Xi 7 T)\ ~* 0- The differences 



18 



will become significant if the initial point is close to the boundary between the 
basins of attraction of different attractors on M., or if the attractor on Ai has 
a positive Lyapunov exponent. However, these are inhospitable circumstances 
to which to expose our technique, and we have shown elsewhere M that for 
a simple atmospheric model the procedure we have described for determining 
the initial slip produces a very good agreement between P(t) and Q(t) over a 
reasonable time-scale for forecasting, even when only the first few terms in the 
power series are kept. 
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A Other invariant manifolds 

In this appendix we generalise the statements we have so far made about systems 
with a centre manifold to other types of invariant manifolds. In the context of 
forming low-dimensional dynamical models this is useful because, for example: 
some systems, when linearised, also have eigenvalues with positive real parts B; 
some decaying modes may not decay fast enough to be reasonably eliminated 
in a given application J23L BQ] . 

A.l Including linearly unstable modes 

First consider the case of a system where some eigenvalues of the linear operator 
have positive real parts: 

x = ,4x + M(x,y,z) xeR m 

y = £y + N(x,y,z) yel" (47) 

z = Cz + L(x,y, z) zGM 1 , 

where the eigenvalues 7^ of C satisfy T > ^R(^k) > 7 > 0, those of A satisfy 
3?(afc) = 0, and those of B satisfy 3?(/3fc) < —ft < 0. The normal form evolution 
equations are of the form 

X = A X + 3-(x)+a{x,V,0 

r, = Bri + B(x,r),0v (48) 

C = cc + c(x,vX)(, 
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where a(0) = a(x, 0, C) = a(%> V, 0) = 0. 

Immediately we identify five non-trivial invariant subspaces of this system 
of equations: 

• the centre space with rj = (, = 0; 

• the stable space with \ — C = 0; 

• the unstable space with \ — V = 0; 

• the centre-unstable space with 77 = 0; 

• the centre-stable subspace with £ = 0. 

In terms of the original variables these five spaces respectively form five invariant 
manifolds: a centre manifold M.\ a stable manifold M-s', an unstable manifold 
M.\j\ a centre-unstable manifold M.cv\ an d a centre-stable manifold, AA.cs- 

In the context of constructing low-dimensional models which capture the 
long-term dynamics of the system, only one of these five invariant manifolds 
is of real interest. Solutions on the centre manifold, the centre-stable manifold 
and the stable manifold arc all unstable to perturbations in the linearly unstable 
variables C, near the origin. Solutions on the unstable manifold, while stable to 
perturbations in 77, may be unstable to perturbations in the centre variables \- 
This leaves the centre-unstable manifold AAcu which is invariant and, at least 
near the origin, attracts all solutions in its neighbourhood at a rate like e"' 3 *. 

On the centre-unstable manifold, the dynamical behaviour is described by 

C = CC + C'(x,0,C)C, x = A X + a(x)- (49) 

Unlike the case without linearly unstable modes, the normal form (|48| ) does 
not in general allow the derivation of initial conditions for a solution Q(t) on 
Mcu that matches a solution P(t) off M-cu- The reason is possible resonances 
between the stable and the unstable modes, which only affect the short-term 
transients (as r\ — > 0) but which cannot be removed from the normal form 
equations (Eq). This is most easily seen in a simple example. 
Consider the normal form equations 

x = v( 2 , v = -2 V , C = C-C 3 - (50) 

These equations have the centre-unstable manifold r\ = which attracts all 
neighbouring solutions as r\ = ?yoexp(— 2t) — > 0. Therefore any neighbouring 
solution P(t) will asymptote exponentially quickly to a solution Q(t) on M.cv\ 
hence for all neighbouring initial conditions P(0) there is a corresponding initial 



condition Q(0) on Mcu- By solving (50) explicitly we find that the solution 
starting from 

X =Xo 1 _ > 2 , v =0, C =Co 
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has the same long-term dynamics as that starting from (xo,?7o,Co)- 

The presence of the logarithm in this expression is novel. It indicates that 
the assumed power-series representation for the normal form transformation is 
not sufficiently general for our purposes. Indeed if we make the nonlinear change 
of variable 

r?C 2 logC 

and leave r] and £ unchanged then (|5Q) becomes 



9 = 0, r) = -2 V , ( = C ~ C 3 ■ 

The new centre variable 9 is not forced by any resonant term, and the projection 
of initial conditions onto the centre-unstable manifold is trivial. 

Let T be a realistic upper bound on 5R(7fc) of C and let — /3 be a realistic 
upper bound on §?(/?&) of B (see text just after (]47|)). Then if T > /3, that is, 
if the most rapidly growing mode grows faster than the slowest decaying mode 
decays, just one of the decaying modes may directly interact with just one of 
the growing modes — the resonance could be as direct as a second-order rnQ 

term. However, if 3? < /3, that is, if the slowest decay occurs more rapidly than 

B/V 
the most rapid growth, then any resonance must be of order T]iQ or higher. 

Thus the normal form projection of initial conditions can be carried out to this 

order of accuracy. Observe that as T — > 0, when the centre-unstable manifold 

becomes a centre manifold, J3/T — > oo, and the order of the resonance moves 

out to infinity. This limit recovers the case of a centre manifold where there is 

no problem in computing an appropriate initial condition Q(0). 

A. 2 Keeping some decaying modes 

In a particular application it may be that some of the linearly stable modes 
of (h_7|) are dynamically significant, for example, if their decay is relatively slow. 
An example is shear dispersion in a channel or pipe, where the slowest modes 
decay on a cross-stream diffusion time-scale that may be comparable with other 
physical processes of interest. We then try to construct an invariant manifold 
Mi that includes not only the critical x-modes together with the unstable modes 
z, but also the slower-decaying of the y-modes. For shear dispersion this was 
done by Watt & Roberts gj. 

After renaming the z variables and some of the y variables of ( j47| ) as x vari- 
ables, we consider a dynamical system given by (fy-0) but where the eigenvalues 
«fe of A and those of B satisfy 

8109*) < -0 < -7 < »(a fc ) < r 0,7 >0. (51) 

Again without loss of generality we take A and B to be diagonal matrices. The 
variables x are termed master variables in this context [Q~3| , while the variables 
y are termed slave variables. 
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As before, we make a nonlinear change of variables. Resonance may now 
generate a term that includes purely master variables as a forcing in the slave 
equations. Thus, even in the normal form, the slaved variables do not necessarily 
evolve quickly to zero. This means that they may evolve non-trivially at large 
times; that is, on the long-term time-scale of the master variables ("short-term" 
then refers to the time-scale of the slave variables' decay). A resonance occurs 
at lowest order when j master modes combine to force one slave mode; using 



the bounds ( 51 ) on the eigenvalues, the lowest order at which the resonance can 
occur is 

js = ph ■ 

This resonance is symptomatic of the non-uniqueness of invariant manifolds. 
The non-uniqueness first appears at this order; as discussed by Roberts [£3| it 
appears as a zero divisor in the construction of the invariant manifold — there 
is little point in computing higher orders in the expansions because the higher 
order differences are smaller than those ignored by the model. 

Resonant interactions between purely master modes are the important inter- 
actions in the long term. However, any interaction with a slave mode results in 
the undesirable feature that during the evolution onto Aii the master dynam- 
ics depend upon the slave dynamics — an effect which we want to ignore. The 
lowest order at which this resonance may appear is if one slave mode interacts 
with j — 1 of the most rapidly growing master modes to force the most rapidly 
decaying master mode. Thus the lowest order for this resonance is 

Such resonance makes it hard to prescribe correct initial conditions on the low- 
dimensional manifold Mi for any given initial condition of the full system. 
However, in many examples, such as (pCf), this reflects the limitations of the 
standard power-series representation. 

Observe that these effects show up at high-order only if the spectral gap 
between the master and the slave modes, [— j3, — 7], is large, both relatively 
(large P/j) and absolutely (large ^jP-), when compared to the natural growth 
rates in the dynamics. Thus such a large spectral gap is desirable in constructing 
low-dimensional dynamical models. 
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